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> '. Introduction 

(N 

Q I Ever since it became possible to perform computer simulations of systems with a few 

hundred particles this method was also used to investigate the structure and dynamics 

t — ■ 

Q>^ I of supercooled liquids and glasses. The absence of long range order makes these systems 

^ ■ a very difficult task for any sort of analytical calculations and since experiments often 

cannot give all the desired information, be it of structural or dynamical nature, com- 
puter simulations offer a very convenient way to access such information. The goal of the 
present article is to review some exemplary results of investigations in which the dynamics 
of supercooled liquids was studied. In this we will focus on recent studies of the proto- 

rS I typical network glassformer Si02, since the lack of space does not permit us to review the 

large body of literature which exists on the general subject of computer simulations of 
supercooled liquids and glasses. However, the interested reader can find a more thorough 
discussion of this subject in the review articles [jl|, 0, ||, Q. 

Before we start to introduce the model we used in our simulations and to present the 
results, we briefiy discuss what the kind of questions are that can be addressed by means 
*To appear in Analysis of Composition and Structure of Glass, and Glass Ceramics Eds.: H. Bach 
and D. Krause (Springer, Berlin, 1998) 
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of computer simulations. Although it is of course in principle possible to make a full 
scale ab initio simulation of any given material, the computational demand restricts such 
types of investigations to systems that are relatively small (less than 100 particles) and 
to short times (on the order of lOps). Such system sizes and time scales are for most 
applications in the field of glasses too small/short and thus one usually limits oneself 
to simulations in which the system is described by a classical force field. However, it is 
essentially impossible to devise classical force fields that reproduce exactly the forces in 
the real material and in practice thus one always uses some approximative potentials. 
Therefore it is clear that it is most improbable that any simulation will reproduce all 
experimentally determined quantities (density, bond-angles, viscosity, density of states, 
etc.). In Ref. 0, e.g., it is nicely demonstrated how various potentials for Si02 can 
give very different answers regarding the absolute numbers of various physical quantities, 
such as the temperature at which the model shows a maximum in the density or on 
the temperature dependence of the diffusion constant, but the bottom line is that most 
reasonable potentials agree at least qualitatively with the experimental findings. Thus, if 
the potential model is sufficiently realistic it can be expected that such a simulation will 
reproduce at least the salient features of the material and that therefore one can use the 
simulation to identify trends (e.g., temperature dependence of the density) or mechanisms 
(e.g., diffusion mechanism) that are present in the real material as well. 

Despite the fact that usually one does not demand a perfect agreement between the 
results of a simulation and real experiments, it is often not a simple task to find a potential 
which gives at least a fair agreement. An example for this situation is B2O3 for which it 
seems to be quite difficult to come up with a good classical potential ^. Thus it is quite 
often the case that the main obstacle to use computer simulations in order to gain insight 
into the properties of supercooled liquids and glasses are not the lack of computational 
resources but the shortage of classical force fields that give a realistic description of the real 
material. Therefore it would be most useful if progress would be made in the development 
of classical force fields. 

Model and Details of the Simulation 

In this section we present the model we used for our Si02 simulation and discuss some of 
the technical details of the simulation. More details can be found in Refs. 0, |^, |^, |TD|. 



A few years ago van Beest et al. introduced a model for Si02 in which the interactions 



between the ions are of a pure two-body type |rT| . Despite the absence of three-body in- 
teractions this model is able to generate the tetrahedral network of silica which is thought 
to give the correct local structure of this system. It was subsequently shown that this 



potential gives a good description of the various crystalline phases of silica [|T^] and of 



the static properties of amorphous Si02 as well [|13| , [14| . Thus it is reasonable to use this 
potential also to investigate the dynamical properties of supercooled silica. The potential 
energy between two ions of type i and j {i,j G {Si,0}), a distance Vij apart, are given by 

' ij ' ij 

The values of the parameters g^, Aij, Bij and Cij can be found in the original publica- 



tion [|TT|. The non-Coulombic part of the potential was truncated and shifted at 5.5A. The 
simulations were done at constant volume and the density of the system was fixed to 2.37 
g/cm^. The system size was 8016 ions, giving a volume of the cubic box of (48.37A)^. 
This size is significantly larger than the ones commonly used in such simulations but 
necessary in order to avoid finite size effects which have been shown to be remarkably 
pronounced in network glassformers ||^. The equations of motion were integrated with a 
time step of 1.6 fs. The temperatures investigated were 6100 K, 5200 K, 4700 K, 4300 K, 
4000 K, 3760 K, 3580 K, 3400 K, 3250 K, 3100 K, 3000 K, 2900 K and 2750 K. At all 
temperatures the length of the runs were longer than the typical relaxation time of the 
system and extended at the lowest temperature over about 20 ns. More details on the 



simulations can be found in Ref. [T^. 



Results 

One of the simplest dynamical quantities that can be obtained from a computer simulation 
is the diffusion constant D of a tagged particle, which can be computed from the long 
time behavior of the mean squared displacement. Care has to be taken, however, that D 
is determined from a time window in which the mean squared displacement has indeed 
reached its asymptotic time dependence, i.e. shows a t^ dependence. Otherwise the 
computed value of D will be too large. 

In Fig. m we show the temperature dependence of Dsi and Dq, the self diffusion 
constant for the silicon and oxygen atoms, respectively. We see that for all temperatures 
investigated Dq is larger than Dsi? which is to be expected since the oxygen atoms are, on 
average, bound to the network by only two covalent bonds, whereas the silicon atoms are 



bound by four bonds. For low temperatures the diffusion constants show the Arrhenius 
temperature dependence well known from experiments. The two straight lines are an 
Arrhenius fit to the low temperature data with activation energies of 4.66eV and 5.18eV for 
the oxygen and silicon atoms, respectively. These values agree well with the experimental 
values determined by Mikkelsen |jl5| and Brebec et al. [|16], see figure, and thus we conclude 



that the used model for Si02 is indeed able to give a quite realistic description of the 
dynamics of supercooled silica. (We note that there is some evidence that for the model 
of van Beest et al. the activation energy for silicon might, at low temperatures, be even 
a bit larger than the value obtained from the presented fit [|l3]-) From the figure we also 
recognize that at higher temperatures the diffusion constants deviate from the Arrhenius 
behavior observed at low temperatures in that they are smaller than expected from an 
extrapolation with the Arrhenius law. So far this behavior has not been observed in real 
experiments, since the temperatures accessed so far in the experiments are still a bit lower 
than the ones for which the non- Arrhenius behavior is observed. But of course we cannot 
exclude the possibility that the observed non- Arrhenius behavior is an artifact of the used 
potential with no counterpart in the real material. On the other hand, the investigation of 
Hemmati et al. has shown that many silica models show such an non- Arrhenius behavior, 
although at different temperatures, thus giving support to the expectation that also real 
silica will show this behavior at sufficiently high temperatures. 

In order to gain some insight why at high temperatures the temperature dependence 
of the diffusion constants seems to be different from the one at low temperatures, we can 
study the dynamics of the particles in more detail. This can be done with the help of the 



self part of the van Hove correlation function which is defined as ||17 



Gi-\r,t) = —(^Y.Kr-m)-f^m)j ,«G{Si,0} . (2) 

Here (.) stands for the thermal average. Thus G^°'\r,t) is the probability that a particle 
of type a has moved within time t a distance r. In Fig. ^ we show 47rr^Gp)(r, t) versus 
r for different times. From these figures we recognize that at high temperatures. Fig. ^, 
the function is a single peaked function with a maximum that is moving to larger values 
of r when time increases. With decreasing temperature the speed with which the peak 
moves to larger distances decreases, but qualitatively the shape of the curves is similar to 
the ones shown in Fig. |^. This behavior is very similar to the one observed for fragile 
glassformers 0, |[ , where at low temperatures the dynamics of the particles is described 



very well by the so-called mode-coupling theory [^, ^, which predicts such a time and 
space dependence of Gs(r, t). If the temperature is lowered further, say below 3500 K, this 
qualitative picture changes, in that for intermediate times Gs{r,t) shows a double peak 
structure, see the curve for t = 4.5 ns in Fig. Qb. This secondary maximum, marked by an 
arrow, is usually associated with the presence of so-called hopping processes, i.e. types of 
movements in which the particle jumps abruptly to a site in the network which is about one 
interparticle distance away from the starting place . These sort of processes are usually 
activated, since the jumping particle has to overcome a barrier formed by its surrounding 
particles (this is in contrast to the motion at higher temperatures, where the dynamics 
of the particles is smoother, i.e. where no jumps occur, see Fig. ||a). Thus the study of 
the van Hove correlation function shows us that the dynamics of the particles changes 
from a smooth movement to a jump-like motion when the temperature is lowered. It is 
this change in the microscopic dynamics that gives rise to the change in the temperature 
dependence of the diffusion constants. 

A further important transport quantity is the shear viscosity rj which can also be 
measured in real experiments. Unfortunately the calculation of t] is not a simple task 
in computer simulations, since it is necessary to average this quantity over quite a few 
(at least 10) a-relaxation times (defined below) before rj is determined with a reasonable 
accuracy. Therefore we were able to measure rj only for temperatures larger than 3000 K. 
The temperature dependence of rj is shown in Fig. ^ in an Arrhenius plot. Despite the 
noise in the data we can recognize that this temperature dependence is not Arrhenius like 
thus is similar to the temperature dependence of the diffusion constant. Also included in 



the figure are the experimental data points of the measurements of Urbain et al. |]20|. We 
see that a reasonable extrapolation of our data to lower temperature will underestimate 
the viscosity at the temperatures at which the experimental values are available, which 
shows that with respect to this quantity the potential of van Beest et al. is perhaps not 
quite reliable. (However, it should also be kept in mind that the experimental data points 
might be subject to an appreciable systematic error.) 

Since experimentally the viscosity is simpler to measure than the diffusion constant, it 
is often common to calculate the latter from the former by means of the so-called Eyring 
equation (which is essentially identical to the Stokes- Einstein equation), i.e. 

D = kBT/7]X. (3) 

Here ks is Boltzmann's constant and the constant A (having the dimension of length) 



is related to the size of an elementary diffusion step and is usually determined from 
measurements at lower temperatures. Having measured both, D and r^, we can now check 
whether the combination ksT /rjD is indeed constant. This quantity is shown in the inset 
of Fig. ^for the silicon as well as the oxygen diffusion constant. From this plot we clearly 
see that in the temperature range investigated the assumption that A is constant is a 
poor approximation. It is, however, interesting that at the lowest temperatures the value 
of A for the oxygen atoms is about 5A, which is not too different from 2.8A, the values 



assumed by Poe et al. |21| at lower temperatures. 

The fact that the atoms of the system form a three dimensional open network implies 
that the bonds between adjacent atoms are broken if an atom diffuses. Therefore it is 
reasonable to look into the time dependence of these bonds. For this we define the quantity 
Pb{t) which is the probability that a bond which is present at time t = is also present at 
time t. (Note that we call a silicon and oxygen atom bonded if their distance is smaller 
than 2.35A, the location of the minimum between the nearest and second nearest neighbor 
peak in the Si-0 radial distribution function. This location is essentially independent of 
temperature.) 

In Fig. ^ we show the time dependence of Pf, for all temperatures investigated. From 
this figure we recognize that at low temperatures the curves do not decay to zero within 
the time span of our simulation. From our investigation of the intermediate scattering 
function, discussed below, we know however, that this correlation function, which mea- 
sures the structural relaxation of the system, does decay within the time span of our 
simulation. Thus we come to the conclusion that for the structural relaxation to occur 
it is not necessary that all Si-0 bonds are broken, but only about 80% of them. This 
number is quite reasonable since it can be argued that an oxygen atom can make a relax- 
ational movement if one or two of its bonds are broken and that a silicon atom can move 
if three of its bonds are broken, i.e. about 75% of the average number of bonds. 

From Fig. ^ we also see that the shape of the curves is essentially independent of 
temperature (a small temperature dependence is observed for times at which P^ is still 



close to its initial value [T^). Thus we come to the conclusion that the mechanism that 



leads to the breaking of the bonds is also essentially independent of temperature. Thus it 
is reasonable to define the mean lifetime rj, of a bond by requiring that Pb{Tb) = e^^. This 
lifetime can now be compared with the typical time for the diffusion. For this we plot in 
Fig. ^ the products DsiTb and -DoTfe as a function of temperature. We recognize from this 
figure that this product is essentially constant for the oxygen atoms, but decreases for the 

6 



silicon atoms. Thus we conclude that the breaking of a bond is related to the diffusive 
movement of the oxygen atom, which is very reasonable, since this atom becomes mobile 
after the breaking of one of its bonds to its two nearest neighbors silicon atoms. For 
silicon atoms the situation is different: Even if one bond to its neighboring oxygen atoms 
is broken it is still tightly bound to the network by the remaining oxygen neighbors and 
thus cannot start to diffuse around. Therefore the quantity Tf, is not related to the diffusion 
constant of the silicon atoms and hence TbDsi is not constant. 

As mentioned above, the structural relaxation of the system can be studied very well 
by means of Fs{q,t), the (incoherent) intermediate scattering function for wave- vector q. 
A further reason why this function is important is that it can be measured in neutron 



and light scattering experiments. It is defined as |T^ 



i^i"ng,t) = ^^Eexp(zg-(r-;(t)-f,(0))^, a G {Si,0} . (4) 

In Fig. I^a we show the time dependence of F^'~'\q, t) for the wavevector q = 1.7A^^ which 
corresponds to the location of the first sharp diffraction peak in the static structure factor. 
We see that at high temperatures the curves decay quickly to zero with a time dependence 
which is approximated well by an exponential function. When the temperature is lowered 
a shoulder starts to form at around 0.4 ps which, upon lowering the temperature further, 
develops into a plateau. The time range during which the correlation function is close 
to the plateau is usually called the /3-relaxation regime, whereas the time range during 
which the correlator starts to fall below the plateau and to zero is called the a-relaxation 
regime. The physical significance of this plateau is that in the time range where it is 
observed the particles rattle around in the cage formed by the neighboring particles and 
only for larger times this cage starts to break up, the particles begin to diffuse and thus 
the correlation function starts to decay to zero. 

From Fig. ^ we also see that at low temperatures the shape of the curves seems, in the 
a-relaxation regime, to be independent of temperature. To check whether this is indeed 
the case we define the a-relaxation time r(g) by requiring that Fs{q, '^{q)) = e~^. Thus, if 
the shape of the curves are in the a-relaxation regime indeed independent of temperature, 
a plot of the curves versus the rescaled time t/T{q) should give a master curve. (If the 
correlation functions show such a scaling behavior one also says that they obey the time- 
temperature superposition principle.) That this is indeed the case is demonstrated in 
Fig. ^, where we plot the same correlation function as shown in Fig. ^ versus t/r(g). 
We see that we find indeed that the curves for the different temperatures fall, for low 



temperatures, reasonably well onto a master curve. Thus we can conclude that the slowing 
down of the relaxation dynamics with decreasing temperature is not due to the fact 
that the decay of the correlation functions at long times becomes slower with decreasing 
temperature but rather that the time until the correlation functions start to deviate 
appreciably from the plateau, i.e. the time for the breaking up of a cage, increases 
strongly with decreasing temperature. 

From Fig. ^one also sees that at times around 0.2-0.5 ps the correlation functions for 
low temperatures show a noticeable dip. Such a feature is not observed in hard-sphere like 



systems, such as Lennard- Jones or soft spheres [g, ^, |2^, which are fragile glassformers, 
and thus seem to be a particularity of strong glassformers. This dip is related to the 
so-called boson peak, a feature which has been observed in various experiments, such as 
inelastic neutron scattering, Raman scattering and more recently also in inelastic x-ray 
scattering. The reason for the presence of this peak is so far a matter of dispute [^ 



smce 



the experiments are not able to provide the information necessary to decide between the 
different theoretical models. Thus this type of question might very well be one which will 
ultimately be decided by means of specific types of computer simulations, since those will 
allow to access the necessary information, and some efforts in this direction have indeed 



been made 24 



The last quantities we discuss in this article are the dispersion relations of the different 
modes. Within a computer simulation it is relatively simple to calculate the wave- and 
frequency dependence of the (collective) dynamic structure factors 

1 f°^ 
Sap{q,uj) = T- Fai3{q,t)exp{iujt)dt a,pe{Si,0} , (5) 

zvr J-oD 

with the coherent intermediate scattering functions 

Fa,{q,t)= l\ M cirexpHg-f)(EE^(r-|r-;(t)-f,(0)|)), (6) 

^^"^^/3 -^ \fc=li=l / 

since all the necessary information to compute this quantity is readily accessible. The 
dynamic structure factors are interesting for at least two reasons: Firstly they can be also 
measured in experiments, although only in a relatively limited range of frequency and 
wave-vectors, and secondly they allow insight into the nature of the excitations present 
in the system, such as the acoustic modes and the optical modes. The fact that we can 
compare the results of Sapiq-, ^) from a simulation with real experimental data gives us the 
possibility to test whether the silica model used is able to reproduce this type of dynamic 
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observable in a satisfactory way. In the case where such a comparison is satisfactory we 
then can thus use the results from the simulation to gain information on Sap{q-,uj) in 
those q and u ranges which are not accessible to real experiments. Such comparisons 
have indeed been made and it was shown that the present models for silica agree quite 
reasonably with experiments [|lOi ^. 



Once the spectra are obtained as a function of frequency the location of the various 
peaks can be read off and thus the dispersion relation of the different modes determined. 
(Note that for the investigation of the g-dependence of the transversalm.odL.es it is necessary 
to compute the transversal current-current correlation function Ct{q,uj), which can also 
be expressed simply by means of the positions and velocities of the particles [ITTH .) 

In Fig. 1^ we show the so obtained dispersion relations for the different modes at 
2900 K 1^, |10[. The data for the Si-Si and the 0-0 correlation are shown in filled 
and open symbols, respectively. Starting at small wavevectors, various branches can be 
observed for each species: Two optical ones in the transverse and longitudinal modes, 
denoted by LOl, L02, TOl and T02, and a longitudinal and transversal acoustic mode 
(LA and TA, respectively). In addition to these modes we find at wavevectors larger than 
0.22A^^ an additional peak in the dynamic structure factor at around 2THz (see Fig. ^). 
From its location at higher values of q, which is always around 2THz and a comparison 
with experimental data [^, 26], we conclude that this excitation is the so-called boson 



peak. The wavevector dependence of the different modes becomes relatively complicated 
when q is large and it has been discussed in greater detail in Refs. ^T^. E.g. we see that 
the longitudinal acoustic mode for the silicon atoms (open triangles pointing downwards) 
show a quasiperiodic g-dependence which reflects the fact that the system has a pseudo- 
Brillouin zone at around 1.4A~^ [^. Presently we do not want to dwell further on the 
discussion of the various curves and just point out one last thing. In Fig. ^ we have also 
included two bold solid lines. These present the linear dispersion behavior of the acoustic 
modes in the viscoelastic limit. Note that the slopes of the lines, which is equal to cl and 
ct, the longitudinal and transversal speed of sound, were not taken as a fit parameter. 



but are the experimental speeds of sound at around 1600 K p6[. From the fact that 
for small wavevectors our data points lie very well on (in the case of the longitudinal 
mode) or at least quite close (in the case of the transversal mode) to these lines with the 
experimental slopes we conclude that our simulations are able to reproduce also this type 
of measurements with good accuracy. 



Conclusions 

The goal of this review was not to make a comprehensive review of all the computer 
simulations of supercooled liquids and glasses. Rather we tried to present some exem- 
plary results in order to demonstrate how computer simulations can contribute to our 
understanding of the dynamics of supercooled liquids. As already mentioned in the Intro- 
duction, the most valuable use of computer simulations is presently not to try to reproduce 
the dynamical features of the material of interest to a very high accuracy. Rather it is 
advisable to use rather simple models which are able to reproduce the essential features 
that one is interested in and to make a careful investigation of the properties of such mod- 
els. As we have demonstrated in this article, even such relatively simple models can give 
a quite realistic description of reality. Due to their simplicity, however, it is possible to 
simulate system sizes and time scales that would not be accessible with more sophisticated 
models. Because with today's computer hardware it is possible to simulate systems that 
are reasonably large, on the order of 50A, for time spans (10-50 ns) that are much longer 
than the microscopic times (1 ps), such simulations can indeed give valuable information 
into the dynamics of glass melts and glasses, such as the lifetime of a bond, the mech- 
anism of diffusion or the dynamics of the particles on the microscopic scale, quantities 
that are experimentally accessible only in a limited way, or not at all. Therefore these 
types of computer simulations are a valuable complement to experimental investigations 
to understand the dynamics of such systems. 
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Figure 1: Temperature dependence of the diffusion constant for the sihcon and oxygen 
atoms (filled squares and open circles, respectively). The straight lines are fits to the 
low temperature data with an Arrhenius law giving the stated activation energies. The 



experimental values for the activation energies (from Refs. |jT5[ and [jT6|) are also given. 
Adapted from Ref. |^. 
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Figure 2: The self part of the van Hove correlation function for the oxygen atoms versus 
distance r for times that are evenly spaced on a logarithmic time axis (see labels), a) 
T = 6100 K, b) T = 2750 K. The arrow marks the location of the secondary peak. From 



Ref. 10 
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Figure 3: Main figure: Temperature dependence of the viscosity as determined from the 
simulation, (filled squares), and from the experiments of Urbain et al, (open circles) 
Ref. [^. Inset: Temperature dependence of ksT/Dr] for the silicon and oxygen diffusion 
constant (circles and squares, respectively). From Ref. [ITU . 
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Figure 4: Time dependence of Pbit), the probability that a Si-0 bond which was present 
at time t = is also present at time t, for all temperatures investigated. The horizontal 
dashed line at e~^ is used to define the lifetime r^ of a bond. From Ref. [llO . 
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Figure 5: Temperature dependence of the product between the lifetime of a Si-0 bond 
and the silicon and oxygen diffusion constant (filled and open symbols, respectively). 
From Ref. [m. 
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Figure 6: a) Time dependence of the incoherent intermediate scattering function Fs{q,t) 
for the oxygen atoms for all temperatures investigated. The value of the wavevector q 
is 1.7A^^ and corresponds to the location of the first sharp diffraction peak in the static 
structure factor. The horizontal dashed line at e~^ is used to define the a-relaxation time 
r. b) Same data as in a) but plotted versus the rescaled time t/r{T). From Ref. |^. 
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Figure 7: Dispersion relations for different modes at T = 2900 K. The filled and open 
symbols correspond to the Si-Si and 0-0 correlation, respectively. The first character in 
the labels (L,T) stands for longitudinal and transversal and the second character (0,A) 
for optical and acoustic. The dispersion relation labeled BP is the one for the boson peak, 
a) All wavevectors considered, b) Dispersion relations at small wavevectors for the 0-0 
correlation. The two straight lines are the expected dispersion relation if the experimental 
velocities of sound are assumed I 



Adapted from Ref. |T3 
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